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A.  Background  and  Motivation 


A.l  Introduction 

Optical  measurement  is  central  to  a  wide  range  of  military  and  commercial  applications.  Leveraging  the  ease  with 
which  optical  sensing  may  now  be  followed  by  sophisticated  digital  computing  a  new  computational  imaging  (Cl) 
paradigm  has  emerged  in  which  algorithmic  degrees  of  freedom  are  co-designed  with  the  more  traditional  physical 
optical  degrees  of  freedom  to  yield  entirely  new  types  of  cameras.  Extended  depth  of  focus,  light-field, 
compressive,  and  tomographic  cameras  are  all  examples  of  this  new  class.  Unfortunately,  progress  in  Cl  has  been 
hampered  by  the  relatively  few  physical  degrees  of  freedom  that  are  traditionally  available  within  optical  design.  In 
this  work  we  have  sought  to  "open  the  aperture"  so  that  a  richer  collection  of  optical  degrees  of  freedom  will  be 
available  for  use  within  the  Cl  paradigm.  Generalized  3D  index  distributions  have  been  our  primary  focus.  Weak 
scattering,  metamaterials,  and  volume  holographic  optical  components  are  treated  as  special  cases.  An  important 
goal  of  this  work  has  been  to  identify  novel  design  and  optimization  strategies  that  enable  us  to  understand  the 
relationship  between  3D  index  distribution  and  optical  function.  A  related  goal  has  been  to  identify  potentially 
useful  non-image  mappings  (i.e.,  optical  systems  that  do  not  attempt  to  replicate  the  object  irradiance  at  the  focal 
plane).  For  example,  an  ideal  digital  imager  may  be  viewed  as  one  in  which  the  optical  aberrations  are  not 
necessarily  absent,  but  rather  one  in  which  the  optical  aberrations  are  computationally  invertible.  The  Cl  paradigm 
facilitates  the  exploitation  of  these  non-traditional  mappings  via  use  of  joint  optimization. 

A.2  Computational  Imaging 

Imaging  is  Mission  Critical  -  Optical  imaging  has  become  an  important  element  within  a  wide  variety  of  military, 
intelligence,  homeland  security,  and  commercial  applications.  It  is  difficult  to  imagine  modern  ISR  without  the 
advanced  EO/IR  imaging  capabilities  that  we  currently  have  at  our  disposal  and  it  would  be  accurate  to  say  that 
imaging  has  become  ubiquitous  throughout  the  DoD.  One  important  driver  of  these  capabilities  has  been  the 
continuous  advances  in  focal  plane  technology.  [1,2]  Sometimes  overlooked  however,  is  the  fact  that  all  of  these 
impressive  performance  improvements  have  come  along  with  another  "game  changing"  feature  -  cameras  are 
now  digital.  A  new  "computational  imaging"  (Cl)  paradigm  has  emerged  in  which  algorithmic  degrees  of  freedom 
may  be  jointly  optimized  along  with  the  more  traditional  physical  optical  degrees  of  freedom  to  yield  entirely  new 
types  of  cameras.  [3,  4]  Recent  demonstrations  that  point  to  the  potential  of  this  paradigm  include  extended  depth 
of  field  imaging  (i.e.,  in  which  a  previously  unavailable  capability  is  enabled  via  joint  optimization),  compressive 
imaging  (i.e.,  through  which  reductions  in  sensor  resources  and  improvements  in  image  quality  are  possible),  and 
3D  imaging  using  the  light  field  camera. [5-7]  Unfortunately,  progress  in  this  domain  has  been  limited  by  the  optical 
degrees  of  freedom  that  are  traditionally  available  to  optical  design:  namely  glass  type,  surface  curvature/form, 
and  working  distances.  One  important  goal  of  the  work  proposed  here  therefore  will  be  to  "open  the  aperture"  on 
optical  design  so  that  a  richer  collection  of  optical  degrees  of  freedom  will  be  available  for  use  within  the  Cl 
paradigm. 

Open  the  Aperture  on  Optical  Design  -  Until  recently  cameras  did  not  employ  focal  planes  or  digital  processing. 
These  early  cameras  were  therefore  required  to  generate  images  in  a  format  that  could  be  directly  consumed  by  a 
human  viewer.  It  has  been  convenient  that  simple  spherical  surfaces  can  be  easily  fabricated  and  will  do  a 
reasonably  good  job  of  forming  images.  For  this  reason  most  of  our  effort  in  the  domain  of  "optical  design"  (OD) 
has  been  directed  toward  fixing  the  imperfections  (i.e.,  aberrations)  that  arise  when  we  attempt  to  employ 
surfaces  between  air  and  various  types  of  glass/plastic  to  form  images.  As  a  result  of  this  historic  convenience  and 
despite  extraordinary  progress  in  FPA  technology,  imaging  optics  and  the  associated  OD  activities  haven't  really 
changed  much  over  the  past  several  hundred  years.  The  OD  community  has  spent  most  of  its  intellectual  energy 
and  design  resources  focused  on  a  very  narrow  region  of  optical  design  space  near  the  "identity  map."  What  this 
means  is  that  the  optical  components  in  a  camera  (i.e.,  whether  a  low  cost  cell  phone  camera  or  a  space  telescope) 
are  generally  designed  to  invert  the  effect  of  propagation.  The  goal  of  these  optical  systems  is  to  create  an  identity 
map  between  the  object  space  and  the  image  space.  The  entire  study  of  aberration  theory  is  an  investigation  into 
the  error  that  we  incur  while  trying  to  realize  the  identity  map  -  essentially  ignoring  all  of  the  other  potentially 
useful  mappings  that  could  be  realized. 
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In  order  to  make  the  previous  assertion  a  bit  more  rigorous  consider  figure  la  in  which  a  conventional  incoherent 
optical  imaging  system  is  depicted.  The  operation  of  this  system  is  well-approximated  by  the  linear  mapping 
g(x,y,z,X)  =  JJJJf(x',y',z',X')  hn(x,y,z,X,x',y',z',X')  dx'dy'dz'dV,  where  f(x,y,z,^)  is  the  object  space  irradiance 
distribution,  g(x,y,z,^)  is  the  image  space  (i.e.,  measurement  space)  irradiance  distribution,  and  hn(x,y,z,X,x',y',z',X') 
is  the  kernel  of  the  linear  operator  that  results  from  a  specific  "optical  configuration"  (i.e.,  specific  choices  for  the 
optical  components  and  distances).  We  may  capture  any  such  optical  configuration  in  the  3D  index  distribution 
n(x,y,z,?i).  Note  that  traditional  OD  will  write  hn(x,y,z,X,x',y',z',^')  =  5(x-x',y-y',z-z',A,-A/)  +  a(x,y,z,^,x',y',z',X')  and 
focus  most  of  its  design  efforts  to  minimize  the  aberration  function  a(x,y,z,X,x',y',z',X')  over  some  range  of  the 
object  space  X,  Y,  Z,  and  A.  For  notational  simplicity  it  will  sometimes  be  convenient  to  consider  the  input  and 
output  irradiance  distributions  as  vectors  in  spaces  determined  by  an  appropriately  selected  basis  for  sampling  the 
underlying  continuous  functions.  [8]  We  will  represent  these  object  and  measurement  vectors  as  f  and  g 
respectively  and  note  that  the  integral  transformation  above  may  be  written  as  a  matrix  transformation  between 
these  discrete  spaces  as  g  =  Hnf.  As  discussed  above,  conventional  optical  design  strives  to  achieve  Hn  =  I,  the 
identity  map,  by  minimizing  the  norm  of  the  aberration  operator  A. 


traditional  design  goal:  h(x,y,z,x',y',z')  =  8(x-x',y-y',z-z') 
actual  kernel:  h(x,y,z,x',y’,z')  =  5(x-x',y-y',z-z')  +  a(x,y,z,x',y',z') 


EDOF  design  goal:  h(x,y,z,x',y',z')  =  h,(x-x',y-y')h4(z-z’)  *  8(x-x',y-y',z-z') 


(a) 


(b) 


Mask  in  Image  Plane 


(c) 


Object 


general  3D  refractive  index 
distribution  n(x,y,z) 


0 


g  =  Hf 


- Y - 

Cl  design  goal:  arbitrary  hfoy^x'.y'.z') 

(d) 


Figure  1:  The  evolution  of  optical  sensing  architectures,  (a)  Traditional  imaging  employs  complex  optical  systems  in  order  to 
approximate  the  identity  operator,  (b)  Extended  depth  of  field  imaging  employs  a  pupil  plane  mask  in  order  to  yield  a  more 
general  Toeplitz  operator  and  provide  additional  functionality  when  combined  with  digital  post-processing,  (c)  Feature-specific 
imaging  offers  a  more  general  class  of  system  operator  but  requires  an  intermediate  image,  (d)  Generalized  3D  index 
distributions  offer  the  potential  to  realize  a  much  larger  class  of  optical  transformations  which  when  jointly  optimized  with 
computational  degrees  of  freedom  will  yield  improved  performance  and  reduced  SWaPC. 

Although  it  is  understandable  from  whence  this  focus  on  the  identity  map  derives,  the  advent  of  digital  imaging 
and  the  potential  for  substantial  post-measurement  processing  has  expanded  the  landscape  of  OD.  We  are 
motivated  to  pursue  designs  in  this  broader  space  because  (a)  optical  components  continue  to  represent  the 
dominant  mass  and  volume  costs  in  most  cameras  and  (b)  isomorphic  image  space  representations  (i.e.,  pixels)  are 
redundant  and  so  make  sub-optimal  use  of  photon,  noise,  and  electrical  power  resources.  We  seek  to  substantially 
improve  both  performance  and  SWaP  by  noting  that  a  FPA  need  not  measure  a  pretty  picture  -  it  only  has  to 
measure  the  information  required  for  the  post-processing  algorithm  to  succeed.  This  insight  together  with  recent 
advances  in  optical  fabrication  (e.g.,  aspheric  surfaces,  non-rotationally  symmetric  surfaces,  3D  GRIN  optics,  and 
metamaterials)  have  made  possible  the  study  of  optical  systems  that  do  not  strive  to  implement  the  identity  map. 
A  natural  evolution  has  recently  taken  place  within  the  Cl  community.  The  first  step  has  been  to  generalize 
traditional  lens  systems  (i.e.,  which  strive  to  achieve  the  identity  mapping  Hn=l)  to  include  various  thin 
transparencies  (i.e.,  masks).  When  a  mask  is  implemented  in  the  pupil  of  an  optical  system  a  more  general  Toeplitz 
or  convolutional  mapping  (Hn=Hi)  is  realized.  The  extended  depth  of  field  (EDOF)  camera  is  an  example  of  this 
simple  form  of  Cl.  [5,  9] 
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The  next  step  in  the  evolution  of  Cl  has  been  to  relax  the  shift-invariant  property  that  emerges  from  these  Toeplitz 
mappings  and  seek  optical  techniques  for  realizing  more  general  shift-variant  mappings.  Ideally,  these  systems 
would  fully  generalize  OD,  enabling  any  arbitrary  Hn  allowed  by  the  laws  of  physics.  We  are  not  there  yet,  but  a 
recent  step  in  this  direction  is  represented  by  the  compressive  or  feature-specific  imaging  (FSI)  system  shown  in 
figure  lc.  [12-15]  Once  again  this  system  employs  a  simple  mask;  however,  by  placing  this  mask  in  an  intermediate 
image  plane,  FSI  is  capable  of  realizing  any  positive-valued  energy-conserving  mapping  in  the  transverse  kernel.  So 
the  system  operator  for  this  system,  Hfsi  derives  from  a  kernel  of  the  form  hn(x,y,z,x',y',z')  =ht(x,y,x',y')ha(z,z') 
where  ht(x,y,x',y')  may  be  any  transverse  kernel  allowed  by  physics  and  ha(z,z')  is  constrained  by  the  imaging  optics 
as  it  would  be  in  a  traditional  design.  In  the  transverse  plane  therefore  we  obtain  g  =  HFsif,  where  the  operator  Hfsi 
realizes  nearly-arbitrary  projections  of  the  object  space.  By  virtue  of  making  nearly-arbitrary  measurements  on  the 
transverse  irradiance  distribution,  FSI  has  been  shown  to  offer  advantages  in  both  image  quality  (e.g.,  in  high  noise 
environments  via  the  multiplexed  advantage)  and  task-specific  performance  (e.g.,  via  removal  of  confounding 
variables  and  noise).  As  in  the  case  of  EDOF  we  find  that  a  willingness  to  admit  solutions  in  which  Hn  *  I  has  opened 
the  aperture  in  OD  space  to  yield  performance  capabilities  that  would  otherwise  not  be  possible. 

The  FSI  system  in  figure  lc  may  be  viewed  as  an  existence  proof  of  both  the  physical  realization  and  the  potential 
performance  advantage  of  a  fully  general  optical  transformation.  The  current  approach  to  FSI  represents  an 
important  step  in  the  evolution  of  Cl  but  suffers  from  two  shortcomings:  (a)  it  is  hardware  inefficient  due  to  the 
use  of  reimaging  optics  and  (b)  it  has  not  been  extended  to  allow  for  programming  of  the  axial  or  spectral  kernels. 
These  two  issues  will  form  a  starting  point  for  the  work  proposed  here.  A  central  question  that  we  will  address  is: 
can  fully  general  optical  maps  be  realized  without  reimaging  (i.e.,  using  a  simple  optical  device  with  non-traditional 
surfaces  and/or  volumetric  index  distribution)?  A  system  concept  for  fully  general  optical  measurements  is  shown 
in  figure  Id  and  employs  a  fully  general  3D  index  distribution.  Note  that  all  of  the  systems  discussed  thus  far  (i.e., 
figures  la-lc)  can  be  represented  as  3D  index  distributions  of  a  very  limited  class;  those  index  distributions  that 
can  be  generated  using  homogeneous  materials  separated  by  slowly-varying  surfaces.  It  is  not  currently  known 
what  additional  benefits  may  be  gained  via  use  of  more  general  3D  index  distributions  and  an  important  goal  of 
the  work  proposed  here  will  be  to  characterize  the  types  of  mappings  that  can  and  cannot  be  realized  via  such  fully 
3D  optical  systems.  The  Cl  paradigm  enables  the  exploitation  of  these  alternate  mappings  via  use  of  joint 
optimization  of  optical  and  algorithmic  variables.  The  overarching  theme  therefore  is  to  investigate  non-image¬ 
forming  optical  systems  in  conjunction  with  sophisticated  post-processing  in  order  to  realize  system-level 
functionality  (e.g.,  either  image  reconstruction  to  form  a  pretty  picture  OR  task-specific  functionality  that  may  not 
require  such  a  conventional  representation  at  any  point  within  the  processing  chain)  that  is  superior  to  that  of  a 
conventional  camera. 

B.  Research  Progress  and  Results 

B.l  -  Task  1:  Develop  a  unified  formalism  for  the  analysis  of  optical  imaging  systems  as 
general  linear  mappings. 

This  task  has  focused  on  developing  design  and  analysis  tools  which  will  quantify  the  connection  between  a  given 
3D  index  distribution  and  the  associated  linear  operator  that  acts  between  object  and  measurement  spaces.  We 
have  pursued  two  main  thrusts:  (a)  Linearized  Design  via  Weak  Scattering  and  (b)  Multi-Linear  Operators  in  Strong 
Scattering. 

B.l.a  Linearized  Design  via  Weak  Scattering 

This  section  describes  a  method  of  using  the  first  Born  approximation  to  model  the  propagation  of  light  through  a 
weakly  scattering  gradient  index  (GRIN)  medium  the  index  of  refraction  of  which  is  invariant  along  the  optical  axis. 
It  is  based  on  the  more  rigorous  split  step  Fourier  method. 

Split-Step  Fourier  Method 

In  the  split-step  Fourier  method,  the  optical  axis  in  the  GRIN  medium  is  divided  into  NM  steps,  of  distance  A z.  For 
each  z  position,  E^m\  an  array  of  length  Nx/  models  the  field  distribution.  The  nth  element  in  the  array  is  the 
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complex  amplitude  of  the  field  at  the  position  x  =  x[n]  =  Ax  (n  -  z  =  mAz,  where  Ax  is  the  spacing 

between  samples  in  the  x  dimension.  The  input  light  field  array,  at  z  =  0,  is  given,  and  each  subsequent  light  field 
+-array  is  computed  based  on  the  previous  one.  For  each  step,  to  simulate  the  light  field's  propagation  over  it,  two 
operations  are  performed  on  the  previous  step's  array: 

E(m)  *  (ANm  o  JPAz)£,(m-1)  ( 1 ) 

The  first  operator,  PAz,  is  a  discrete  implementation  of  the  Fresnel  diffraction  integral  over  the  distance  A z: 

gink0Az  (2) 

PAz[m,n\  =  — - ein[x[m]-x[n])2/(AAz^Xm 

el*\[xKz 

where  n  is  the  nominal  value  of  the  index  of  refraction;  A  is  the  wavelength  of  the  light;  and  k0  =  2n/A.  The 
second  operator,  ANm,  accounts  for  the  scattering.  It  is  the  diagonal  matrix  such  that: 

AN m[n,n]  =  elk° AzSnm[n]  (3  j 

where  the  index  of  refraction  at  the  coordinates  (z,x)  =  ( mAz,nAx )  is  n  4-  Snm[n\.  Consequently,  the  field 
distribution  at  the  output  of  the  GRIN  medium  can  be  expressed  as  a  function  of  the  input  light  field  to  the  GRIN 
medium  as  follows: 

£0VM)  =  ANjVm  0  pAz  0  ANjVm_i  0  pAz  o  ...  o  AN,  o  Piz£®  (  4 ) 


Weak  Scattering  Approximation 

If  the  index  of  refraction  in  the  GRIN  is  invariant  in  the  z  dimension,  then: 

AN,  =  -  =  ANWm  =  AN  (5) 

=>  E(Nm)  =  (AN  o  PAz)nmE(0)  (6) 

Let  M  =  (AN  °  PAz)  -  I,  where  /  is  the  identity  matrix.  Then: 

Nm  (7) 

E° Vm)  =  (I  +  M)Wm£(0)  =  V 

m= 0 

The  weakly  scattering  approximation  might  be  made  if  ||M||  «  1.  In  this  approximation,  only  the  first  two  terms  of 
the  summation  are  used  to  compute  the  field  at  the  output  of  the  GRIN: 

£(wm)  «  (/  +  NmM)Ew  (8) 


Example:  GRIN  Lens 

The  example  of  a  GRIN  lens  is  used  to  illustrate  the  validity  and  limitations  of  this  weakly  scattering  approximation. 
The  GRIN  lens  has  an  index  distribution  that  is  a  function  of  x: 


n(x)  =n  — 


(9) 


2fLz 


where  /  is  the  back  focal  distance,  and  Lz  is  the  thickness  of  the  GRIN  lens  along  the  z  dimension.  Flence 

x[n]2  (10) 


Upon  substitution  of  (10)  into  (3), 


Snm[n]  =  —  Vm 

Z/Lz 


x[n ]2 

AN  [n,  n]  =  e  °2fNM 


(11) 


If  the  lens  is  very  thin,  and  then  the  split-step  Fourier  method  calculates  that  the  field  at  the  rear  face  of 

the  GRIN  lens  is  =  AN NmE^°\  If  E^  represents  a  normally  incident  plane  wave,  i.e. 

E (°)  [n]  =  1  V  n  ( 12 ) 

then  the  field  emerging  from  the  lens  is  a  spherical  wave  that  will  converge  a  distance  /  behind  the  lens: 

(13  ) 

E(NM)[n]  =  e  0  2/  ' 

(If  the  lens  is  thicker,  the  back  focal  length  can  be  found  using  Gaussian  reduction,  treating  each  of  the  cross 
sections  as  a  thin  lens  with  focal  length  NMf). 

Meanwhile,  the  weakly  scattering  approximation  calculates  the  field  at  the  rear  face  of  the  GRIN  lens  to  be: 

E° Vm)  «  (/  +  NmM)Ew  =  1 ?(0)  +  NmANE(0)  ( 14 ) 

Therefore,  if  the  field  incident  on  the  front  face  of  the  lens  is  the  normally-incident  plane  wave  described  by  (12), 
the  weakly  scattering  approximation  computes: 
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(15) 


x[n]2 

=>  f  (%)[„]  =  l  +  wMe  Ifco2/wM  =  i  +  wMe-;°[n] 
where  0[n]  =  fc0  .  The  greatest  neglected  term  of  the  expansion  in  (7)  is  AN2#®,  the  nth  component 

2 fNM  \  Z  / 

of  which,  given  (12),  is  So,  for  instance,  in  the  case  where  NM  =  2,  in  order  for  this  weakly 

scattering  approximation  to  yield  accurate  results  for  the  field  at  the  rear  of  the  GRIN  lens: 

\e-i29[n]\  «  !  <=>  2 0[n]  <  7  »  |*[n]|  < 

4  2 

For  example,  Figure  la  and  Figure  lb  plot  fields  computed  with  the  split-step  Fourier  method  and  with  the  weak 
scattering  approximation  for  the  case  where  NM  =  2,  X  =  .5 [im,  f  =  3m.  As  predicted  by  (16),  the  Born 
approximation  (weakly  scattering)  calculations  agree  with  the  split-step  Fourier  calculations  when  |x|  <  .6 mm. 
Figure  2  plots  the  amplitudes  of  the  two  computed  fields  of  Figure  1  after  an  operator  has  been  applied  to  them  to 
implement  the  Fresnel  diffraction  integral  to  propagate  them  3  m  to  the  focal  plane:  The  weak  scattering 
approximation  is  "noisier,"  due  to  the  loss  of  the  higher  order  term. 


Amplitude 


Figure  1A:  \ E^\  computations  with  split-step  Fourier 
method  and  with  weak  scattering  approximation 


Phase,  z  =  0  mm 


Figure  lb:  computations  with  split-step  Fourier 

method  and  with  weak  scattering  approximation 


Amplitude:  z  =  3000  mm 


Figure  2:  Fields  at  Rear  of  GRIN  Lens  Propagated  to  the  Focal  Plane 
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B.l.b  Multi-Linear  Operators  in  Strong  Scattering 

This  sub-task  treats  a  gradient  index  (GRIN)  medium  in  the  shape  of  a  rectangular  prism.  The  objective  of  this  work 
is  to  find  the  index  distribution  of  the  GRIN  medium  such  that  a  certain  set  of  input  light  fields,  incident  on  the 
front  face,  emerge  from  the  rear  face  matching  as  closely  as  possible  to  a  set  of  corresponding  desired  output  light 
fields.  In  order  to  find  this  index  distribution,  a  gradient  descent  algorithm  is  used.  For  each  iteration  of  the 
gradient  descent  algorithm,  the  split-step  Fourier  method  is  used  to  model  the  propagation  of  the  specified  input 
fields  through  the  medium  and  compute  the  distribution  of  the  light  fields  at  the  rear  face  of  the  medium.  The 
output  light  field  distributions  are  compared  with  the  desired  output  light  fields,  and  a  cost  metric  is  calculated  as 
a  measure  of  how  closely  the  realized  outputs  are  to  the  desired  ones.  Based  on  the  gradient  of  this  cost  metric 
with  respect  to  variables  that  are  directly  related  to  the  index  of  refraction  distribution,  the  index  of  refraction  of 
the  medium  is  iteratively  updated. 


The  Split  Step  Fourier  Method 

The  split-step  Fourier  method,  illustrated  in  Figure  1,  is  used  to  simulate  the  propagation  of  the  light  fields  through 
the  GRIN  material.  This  method  of  "beam  propagation"  can  accurately,  model  strong  volumetric  scattering.  An 
Ny  x  Nx  matrix  Uin  represents  discrete  spatial  samples  of  an  input  field,  in  the  plane  z  =  0.  For  the  analysis,  this 
matrix  is  up-sampled  initially  by  a  factor  of  N0,  and  yielding  an  Nyo  x  Nxo  matrix  U0~,  where  Nxo  =  N0NX,  and 
Nyo  =  N0Ny.  U0-  is  also  scaled  so  that  it  has  the  same  norm  as  Uin.  The  optical  axis  through  the  GRIN  material  is 
divided  into  NM  small  steps  of  distance  d .  For  each  step,  two  operations  are  performed  on  this  field  matrix,  to 
simulate  the  field's  propagation  over  the  step.  First,  to  account  for  the  effect  of  deviations  in  the  index  from  the 
nominal  value,  each  value  of  the  field  matrix  is  multiplied  by  a  unit-amplitude  complex  number. 

Um+  \Py>  ^x\  ^m\jly>7lx\  ^  d rri~\jAy,Tlx\  (17) 

where  Um-  is  the  field  matrix  in  the  plane  z  =  md;  Um+  is  the  updated  field  matrix;  and  Mm  is  an  Nyo  x  Nxo 
"mask"  matrix  of  unit-amplitude  complex  numbers: 

Mm[ny,nx\  =  exp(i9m[ny,nx])  (is) 

The  phase  delays  0m  model  the  effect  of  the  deviations  in  the  index  of  refraction  from  the  nominal  value  in  the 
plane  z  =  md.  The  second  operation  performed  to  simulate  the  field's  propagation  over  a  step  is  a  discrete 
approximation  of  the  Fresnel  diffraction  integral. 

U(m+ 1)“  =  (3PUm+Q  (19) 

where  Um+  is  the  field  matrix  after  each  element  has  been  multiplied  by  the  corresponding  element  in  the  "mask" 
matrix;  and  U^m+1y  is  the  new  field  matrix  describing  the  field  in  the  x-y  plane  for  z  =  (m  +  1  )d.  The  scalar  /?  is: 

eikd 

P=—-jAxAy  (20) 

iAd 


where  Ax  is  the  x-spacing  of  the  oversampled  field  matrices;  Ay  is  the  y-spacing  of  the  oversampled  field  matrices; 


A  is  the  wavelength;  and  k  is  the  wavenumber,  — ,  where  A  is  the  wavelength  corresponding  to  the  nominal  index 

A 

of  refraction.  P  is  the  Nyo  x  Nyo  matrix  such  that: 

P[n'y,ny]  (21) 

Q  is  the  Nxo  x  Nxo  matrix  such  that: 


Q[n'x,nx] 


e'2dAx2(n*_n*)2 


(22) 


This  operation  is  a  discrete  approximation  of  the  Fresnel  diffraction  integral.  Upon  substitution  of  (4),  (5),  and  (6) 
into  (3): 


Ut 


(m+l)' 


[n'y,n'x] 


Akd 


iAd 


AxAy  ^ 


,iTdAy2(n'y-nyy 


i 


Um+  \jly 


n x\e  id 


i^x2(n'x-nx)2 


(23) 
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Rearrangement  of  (7)  yields: 


Ny  Nx 


!W r[»yXl=7jj£  Z 

ny=l  nx=l 

Therefore,  f/(m+1)-  [n^n*]  is  a  discrete  approximation  of  the  Fresnel  diffraction  integral: 

oo 

U(m+i)-(n'xAx,n'yAy)  =  6—  JJ  Um+(f,T])el2 +(n'y*y-n)  \d^drj 


(24) 


(25) 
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Figure  3:  Split  Step  Fourier  Method 

The  goal  of  this  optimization  is  to  find  the  phases  of  the  "mask"  matrices,  Om[ny,nx]  for  1  <  m  <  Nm,  1  <  ny  < 
Ny,  1  <  nx  <  Nx  that  minimize  the  difference  between  the  field  matrices  that  the  system  outputs  when  the 
specified  input  field  matrices  are  sent  through  it,  and  the  desired  output  field  matrices. 

The  Gradient  Descent  Algorithm 

Each  iteration  of  the  gradient  descent  algorithm,  the  phases  of  the  masks  are  updated: 

drn +1)[%. nx]  =  CK  nx]  -  r]Vdmcost[ny, nx]  (26) 

where  d^[ny,nx]  is  the  element  in  row  ny  and  column  nx  of  9m  on  the  kth  iteration;  d^+1->[ny,nx]  is  the 
element  in  row  ny  and  column  nx  of  6m  on  the  k  +  1st  iteration;  Vemcost[ny,nx]  is  the  partial  derivative  of  the 
cost  metric  with  respect  to  Om[ny,nx];  and  77  is  a  scale  factor. 

The  cost  function  used  by  the  gradient  descent  is  a  weighted  combination  of  the  squared  error  between  the 
desired  and  realized  outputs,  and  a  function  of  the  correlation  between  the  desired  and  realized  outputs.  Let  U^t 
denote  the  output  field  matrix  realized  by  applying  the  split  step  Fourier  propagation  method  to  the  nth  specified 
input  field  matrix,  and  let  U^tdes  be  the  desired  output.  Then  the  cost  associated  with  the  nth  input/output  pair 
is: 

cost^  =  Atx  trace  (^U^t4es  -  U™)\u™4es  -  U^t)j  +  A2p  (27, 

where 

|trqCe(^)ttf/0(")Mes)|2 

trace(U^l4Ju(^t4es)  x  traced U™) 

and  where  X1  and  X2  are  weights.  The  total  cost  for  one  iteration  is  the  cost  averaged  over  the  NI0  specified 
input/output  pairs. 
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(29) 


cost  = 


£^7°  COSt (n) 

n7o 


Example:  Simple  Lens  Optimization 

In  this  example,  the  desired  input/output  field  matrix  pairs  "point-source"  input  field  distributions  to  "point- 
source"  output  field  distributions,  as  a  lens  would.  Nx  =  Ny  =  6.  The  other  parameters  of  the  optimization  are 
listed  in  the  table  below. 


Table  1:  Parameter  Values  for  Simple  Lens  Optimization 


Parameter 

Value 

Nx 

6 

Ny 

6 

No 

20 

A 

.  5/1771 

d 

1 771771 

Nm 

6 

h 

.1 

a2 

1 

Nit  (#  of  iterations) 

1000 

V 

1 

There  are  36  specified  input/output  pairs.  The  nth  input  field  matrix  is  zero  everywhere,  except  that  it  is  1  in  the 
row  r  and  column  c  where  r  =  mod(n  —  1, 6)  4-  1,  and  c  =  floor  (““)  +  1-  The  nth  corresponding  desired 
output  field  matrix  is  zero  everywhere  except  in  the  row  7  —  r  and  column  7  —  c.  For  instance,  the  13th 
input/output  field  matrix  pairs  are  illustrated  in  Figure  2. 


Input  Field  Matrix  Desired  Output  Field  Matrix 


123456  123456 


Figure  4:  Simple  Lens  Input/Output  Field  Pair 

These  desired  outputs  correspond,  roughly  to  the  output  fields  that  a  simple  lens  would  produce  from  the  given 
input  fields. 

The  optimization  is  run  5  different  times.  Each  time,  a  different  set  of  the  "mask"  matrices  are  allowed  to  be  varied 
by  the  gradient  descent  algorithm,  while  the  rest  of  the  "mask"  matrices  are  fixed  at  all  l's,  corresponding  to  no 
deviation  in  the  index  of  refraction  from  the  nominal  value.  The  different  sets  of  mask  matrices  that  are  allowed  to 
vary  are  listed  in  Table  2. 
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Table  2:  Mask  Matrices  That  Vary 


Trial 

"Mask"  Matrices  That  Are  Varied 

1 

m3 

2 

m2,m4 

3 

4 

m1,m2,m4,m5 

5 

mvm2,m3,m4,m5 

This  experiment,  furthermore,  is  run  8  different  times,  each  time  initializing  the  phases  of  the  "mask"  matrices 
differently.  In  7  of  the  8  cases,  the  phases  of  the  "mask"  matrices  are  initialized  randomly,  with  a  uniform 
distribution  between  0  and  2n,  and  in  the  8th  case,  the  "vacuum-initialization"  case,  all  phases  of  the  "mask" 
matrices  are  initialized  as  0's. 

As  the  trial  index  increases,  the  number  of  free  "mask"  matrix  parameters  increases.  As  this  number  of  degrees  of 
freedom  increases,  it  is  expected  that  the  final  cost  metric  and  the  RMSE  between  the  realized  field  matrices  and 
desired  field  matrices  decreases. 


Figure  4  plots  the  cost  metric  as  a  function  of  the  number  of  "mask"  matrices  being  varied.  Specifically,  the  blue 
curve  represents  the  average  cost  of  the  7  randomly-initialized  trials;  the  orange  and  yellow  curves  represent  the 
average  cost  plus  or  minus  a  standard  deviation.  The  purple  curve  represents  the  "vacuum  initialization"  case. 
While  the  cost  in  the  "vacuum  initialization"  case  decreases  monotonically  with  the  number  of  varying  "mask" 
matrices,  as  expected,  the  cost  increases  on  average  for  the  randomly-initialized  cases. 


Cost  v.  #  Masks 


-  /z 

—A - 

/Z  C 7 

- 

fi  a 

— X- 

vacuum-initialized 

#  Masks 

Figure  5:  Cost  v.  #  Masks 


This  trend  is  verified  by  Figure  4a-c,  which  plots  U^ut  for  the  vacuum  initialization  case  with  one  variable  "mask" 
matrix,  for  the  vacuum  initialization  case  with  5  variable  "mask"  matrices,  for  the  random  initialization  case  with 
one  variable  "mask"  matrix,  and  for  the  random  initialization  case  with  5  variable  "mask"  matrices.  The  case  of 
random  initialization  case  with  5  variable  "mask"  matrices  appears  worse  than  the  others. 


9 

Approved  for  public  release:  distribution  is  unlimited. 
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( 13l 

Figure  B.l.b.6a:  JJout  Vacuum-Initialized 1  mask 


Figure  B.l.b.4b:  Vacuum-Initialized ,  5  masks 


/VtfualOufcut 


20  40  so  m  iaoi20 


20  40  60  SO  100120 


Figure  B.l.b.4c:  U^j)  Randomly-Initialized,  1  mask 


Figure  B.l.b.4d:  Randomly-Initialized,  5  masks 


Moreover,  the  RMSE  increases  with  the  number  of  free  "mask"  matrices.  Figure  5  plots  the  mean  RMSE,  as  a 
percentage  of  the  dynamic  range,  versus  the  number  of  variable  "mask"  matrices.  Specifically,  the  blue  curve 
represents  the  average  RMSE  of  the  7  randomly-initialized  trials;  the  orange  and  yellow  curves  represent  the 
average  RMSE  plus  or  minus  a  standard  deviation.  The  purple  curve  represents  the  "vacuum  initialization"  case.  It 
is  hypothesized  that  this  increase  in  cost  and  RMSE,  for  the  randomly  initialized  cases,  is  due  to  an  increase  in  local 
minima  in  the  cost  function  as  the  number  of  "mask"  matrices  increases. 


RMSE  v.  #  Masks 


Figure  7:  RMSE  v.  #  Masks 
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B.2.  -  Task  2:  Investigate  volume  holography  (VH)  as  a  tool  for  the  realization  of  generalized 
non-identity  optical  mappings. 

This  task  has  focused  on  a  quantitative  assessment  of  the  capabilities  and  limitations  of  VH  as  an  optical  element 
for  the  realization  of  arbitrary  mappings.  Kickoff  meeting  feedback  from  DARPA  suggested  that  we  reduce  our 
attention  to  this  task  and  as  a  result,  we  have  only  examined  the  basic  ID  VH  implementation  of  feature-specific 
(i.e.,  compressive)  imaging  using  a  few  well-known  projection  operators. 

Volume  holography  (VH)  is  used  to  develop  a  quantitative  assessment  of  the  capabilities  and  limitations  when  used 
as  an  optical  element  for  the  realization  of  arbitrary  linear  mappings.  It  will  be  almost  trivial  to  capture  the 
behavior  (i.e.,  as  a  linear  operator)  of  each  volume  grating  in  terms  of  its  action  on  the  input  field,  assuming  "ideal" 
cases  of  unit  efficiency,  perfect  Bragg  selectivity  with  no  crosstalk,  no  chromatic  effect,  and  no  degeneracy.  So  our 
task  is  to  capture  these  various  departures  from  ideality  in  the  linear  operator  framework.  In  the  computational 
framework,  the  grating  dimensions  or  thicknesses  are  included  as  a  free  parameter  so  that  the  VH  model  scales 
smoothly  from  thin  (i.e.,  non-Bragg)  to  thick  (i.e.,  Bragg)  behavior.  The  average  refractive  index  n,  modulation 
amplitude  An,  modulation  vector  Kg  of  gratings  are  also  included  as  design  parameters  for  the  optics.  A  volume 
optics  is  suggested  as  a  device  to  support  a  generalized  class  of  computational  imaging.  In  order  to  reveal  the 
connection  between  the  refractive  index  (Rl)  distribution  n(x,y,z)  and  the  system  operator  H,  the  VH  is  used  as  a 
functional  tool  to  analyze/characterize  the  imaging  performance  of  Rl  distribution.  The  Rl  distribution  is  linearly 
decomposed  by  volume  holographic  modulation 


M  -M 

n(x,y,z)=im  ^  T , 

q=\ 

where  N  is  the  coefficient  of  the  VH  modulation  and  q  is  the  number  of  grating.  Note  the  sinusoidal  expansion  of  Rl 
distribution  can  apply  the  diffraction  theory  to  each  component  and  characterize  the  diffraction  properties  of  the 
optics.  Since  any  smooth  function  can  be  written  as  a  superposition  of  sinusoids,  so  both  of  these  pictures  are 
correct  but  we  need  to  be  careful  which  picture  is  physically  relevant.  For  example,  when  the  gratings  have  high 
spatial  frequency  then  diffraction  is  the  dominant  effect.  In  the  other  case,  ray  tracing  is  most  relevant  to 
understand  the  functionality  of  the  optics.  The  VH  grating  produces  a  small  phase  modulation  to  the  dielectric 
constant  s.  The  magnitude  of  the  spatial  modulation  ei  is  much  smaller  than  the  average  magnitude  s  and  the 
modulation  direction  and  period  is  determined  by  the  grating  vector  Kg.  The  total  spatial  modulation  8t  is 
expressed  by  summing  up  the  phase  distribution  of  the  multiple  gratings  with  an  index  k. 

M 

£t  =  £  ■ +  £1  k  ■  cos  (tfg*  •  r) , 
k 

where  r  indicates  the  spatial  position  vector.  Grating  modulation  is  defined  by  the  finite  extent  with  modulations  in 
the  3D  space  coordinates  (x,y,z) 


n(x,y\%)  =  AnRect(^)  Rect(y^}Rect{^-)e 
LX  Ly  Lz 

where  An  is  the  modulation  amplitude  in  the  medium  and  the  exponential  term  is  the  grating  modulation  with  the 
grating  vector  Kg.  To  simplify  the  efficiency  model  of  volume  hologram,  we  start  with  the  ID  grating  with  infinite 
extents  in  the  other  dimensions.  The  ID  grating  uncertainty  is  obtained  by  Fourier  transforming  the  Rl  modulation. 
Assuming  the  infinite  extents  in  the  other  dimensions,  we  ignore  the  Bragg  condition.  Because  of  the  finite  extent 
in  the  dimension,  the  uncertainty  is  side-lobed  by  the  Sine  function  which  determines  the  diffraction  distribution  in 
the  imaging  range. 

n(AK)  =AnLSinc(*^), 
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The  diffraction  distribution  of  volume  grating  is  characterized  by  the  Bragg  condition  and  the  grating  uncertainty  in 
the  FOV.  The  grating  vector  is  determined  by  the  incident  field  vector  Ki  and  the  diffraction  field  vector  Kd  for  a 
single  grating  case  Kg  =  Kd  -  Ki.  In  ideal  thick  holograms,  the  Bragg  condition  is  only  possibility  to  diffract  the  light 
through  the  hologram  medium.  In  practice,  we  have  a  range  of  incident  fields  with  the  vector  Ki  and  the  grating 
uncertainty  relaxes  the  Bragg  condition.  The  diffraction  field  vector  Kd  is  determined  by  the  grating  vector  Kg  and 
the  incident  field  vector  Ki  as  shown  in  Fig.  2. 


(b)  grating  reconstruction 
Fig.  2:  Bragg  mismatch  in  k-space 


At  the  incident  and  diffraction  fields,  the  Bragg  mismatch  AK  is  estimated  by  the  k-space  distance  from  the  point  G' 
to  the  Bragg  circle.  The  Bragg  mismatch  directly  affects  the  side-lobes  of  the  Sinc()  function. 


MC  =  \R0-Rd\, 

The  diffraction  field  with  the  vector  K'g  is  determined  by  the  Bragg  selectivity  in  the  2D  k-space  in  Fig.  3.  Only 
diffraction  which  satisfies  the  Bragg  selectivity  in  the  x  and  z  analysis  practically  occurs  in  the  VH  imaging.  The 
grating  uncertainty  AK  determines  the  efficiency  bands  in  the  ID  imaging  and  the  grating  uncertainty  is 
parameterized  by  the  medium  thickness.  Two  thicknesses  of  1mm  and  lOjum  were  compared  to  show  the  Bragg 
selection  and  the  thinner  grating  showed  the  broader  efficiency  band  with  weaker  diffraction  efficiency. 


Fig.  3:  Bragg  selectivity  in  2D  k-space  analysis:  the  grating  uncertainty  AA'  determines  the 
efficiency  bands  in  the  ID  imaging  and  the  grating  uncertainty  is  parameterized  by  the  medium 
thickness.  (1  mm 2  vs.  10  {im2) 

ID  VH  imaging  is  designed  to  show  a  tangible  example  for  feature  specific  imaging  in  Fig.  4. In  the  operational 
setup,  we  think  about  object  space  being  mapped  to  measurement  space.  Considering  object  space  at  infinity  the 
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VH  analysis  clear  to  explain  because  the  natural  input  modes  will  be  plane  waves  and  these  of  course  will  naturally 
arise  from  object  points  at  infinity. 


0* 


^±20* 


0‘ 


(a)  Imaging  setup 


(b)  k-space 


Fig.  4:  ID  VH  imaging  for  FSI:  90  deg  diffraction  is  the  nominal  mapping  through  the  volume 
medium. 


In  the  measurement  space,  we  need  a  lens-like  element  somewhere  in  the  volume  optics  (i.e.,  as  a  separate  GRIN 
device  or  as  a  quadratic  phase  distortion  on  your  gratings)  to  convert  diffracted  plane  waves  back  to  points  at 
some  finite/small  distance  from  the  GRIN.  Eventually  we  will  be  interested  in  incoherent  illumination,  but  for  now 
we  assume  coherent  but  not  necessarily  monochromatic  illumination.  After  the  focusing  lens,  the  grating 
uncertainty  becomes  the  PSF  for  the  system.  The  linear  operator  is  defined  by  the  shift -variant  PSFs,  which  are 
characterized  by  Bragg  selectivity  and  Grating  uncertainty.  The  grating  uncertainty  is  separately  expressed  in  the  x 
and  z  axes. 


where  the  Bragg  mismatches  are  denoted  by  AKX  =  Kx  -  Kgx  and  AKZ  =  Kz  -  Kgz  in  the  axes. 
The  k-space  distances  are  calculated  by  trigonometry  in  the  k-space. 


A Kz  =  (i  (cos  Bj  +  cos  Go  —  cos  Gj  —  cos  Gj), 
AKX  =  /I (—  sin  Gi  —  sin Gq  +  sin Gi  -b  sin  Gj ) , 


where  p  is  the  average  propagation  constant  and  0  is  the  angle  for  the  corresponding  grating  vectors.  Diffraction 
efficiency  is  directly  related  to  the  uncertainty  distance  in  the  k-space  and  ID  diffraction  efficiency  is  obtained  by 
multiplying  the  individual  diffraction  efficiencies. 


mo  =  nA^d)  ■  T]z(Ri,Rd), 


The  Rl  distribution  n  is  modeled  by  the  linear  combination  of  VH  gratings  thus  the  system  operator  H  is  directly 
defined  by  the  characteristics  of  VH  gratings.  Once  we  find  an  ideal  system  matrix  Hi  from  FSI  imaging,  the  system 
performance  is  characterized  using  the  VH  analysis.  The  system  matrix  Ht  is  estimated  by  summing  the  single 
grating  matrices  Hi  with  proper  weights  wi.  Since  the  system  matrix  is  optically  encoded  by  the  VH  gratings,  the 
system  matrix  is  possible  to  be  made  by  a  linear  combination  of  grating  matrices  Hi.  The  corresponding  weights  wi 
is  obtained  by  projecting  the  ideal  system  matrix  Hi  in  the  space  of  Hi.  Due  to  the  characteristics  of  individual 
gratings  in  the  system  matrix,  the  matrix  Ht  is  not  exactly  equal  to  the  ideal  matrix  Hi ,  which  causes  the 
performance  degradation  in  FSI  applications. 
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M 


H[  +  error. 


For  the  system  operator  of  FSI,  we  examine  Fladamard,  PCA,  and  random  Bernoulli  operators.  The  Fladamard  code, 
composed  of  ‘-V  and  T  transmittance,  is  an  orthogonal  matrix  balancing  the  two  entry  values.  Using  the 
multiplexing  benefit,  the  MSE  is  reduced  by  a  factor  of  4x.  The  PCA  code  is  optimum  in  the  sense  that  they  provide 
dimensionality  reduction  with  the  smallest  possible  MSE.  The  column  vectors  are  principal  components  trained 
from  5  differently  categorized  images.  A  random  Gaussian  code,  generated  by  0  mean  and  1  std,  is  QR  factorized 
to  make  the  orthogonality  in  columns.  A  random  Bernoulli  code,  composed  of  ‘-V  and  T  transmittance,  is 
designed  by  the  Bernoulli  distribution  with  50%  probability.  Thus  the  code  has  the  half  of  opening  and  the  half  of 
closing.  In  compressive  sensing  the  system  matrix  Ht  will  not  be  full  rank  so  the  image  reconstruction  will  require 
some  prior  and/or  regularization.  The  simple  matrix  inversion  H-i  is  impossible,  which  should  reduce  to  Ht  in  the 
noise-free  case.  The  linear  measurement  model  is  vectorized  to  incorporate  with  a  sparse  optimization  and  we  use 
the  vectorized  data  of  the  measurement  g  and  the  object  f  with  the  designed  system  matrix  FI.  The  LMMSE  is 
adopted  by  an  optimal  linear  reconstruction  operator  depending  on  no  free  parameters  and  the  linear  operator 
serves  as  a  well-defined  baseline  for  the  study.  In  the  compressive  measurement  of  g,  the  operator  W  is  defined  by 
the  object  autocorrelation  matrix  R  and  the  noise  covariance  matrix  D.  The  reconstruction  f*  is  obtained  by 
multiplying  the  measurement  with  the  operator, 

g  [  =  fff+n, 

W  =  RH$ -(HRH$+D)-\ 

where  singular  value  decomposition  is  used  to  invert  the  matrix.  The  object  autocorrelation  can  be  obtained  from 
training  data  or  from  one  of  the  many  popular  models  for  natural  images  (i.e.,  the  power  spectrum  of  natural 
scenes).  The  noise  covariance  does  not  need  to  be  trained  it  is  diagonal  (if  iid  noise).  If  the  noise  is  AWGN  then  it 
will  be  proportional  to  the  identity;  whereas,  if  the  noise  is  shot,  then  the  diagonal  will  be  proportional  to  the 
object.  To  decode  the  compressed  measurement  in  g,  the  sparse  optimization  is  used  to  alleviate  the 
underdetermined  problem.  The  system  matrix  FI  represents  a  highly  underdetermined  system,  which  makes  the 
inverse  problem  hard.  Flowever,  an  optimization  procedure  that  imposes  sparsity  constraints  on  f  can  return  the 
most  plausible  one  consistent  with  the  given  g.  The  sparsity  constraints  of  compressive  sensing  [12,  13]  turn  the 
inverse  problem  into  one  that  can  be  solved  using,  for  instance,  11-minimization.  The  sparsity  based  optimization 
can  be  expressed  by 

r  =  argnuni||g  — H  .  f||| +  T||f ||if 

where  f*  indicates  the  estimate  of  data  from  the  optimization  problem,  and  x  is  a  regularization  parameter  for  the 
total  variation  term.  In  the  computational  processing,  TV  based  TwIST  algorithm  was  used  for  the  inversion 
problem. 

Five  test  signals  were  randomly  chosen  to  represent  the  smooth  objects  in  the  ID  imaging  scenario.  Sensing  rates 
of  [5,  25,  45,  65,  85,  100]  %  were  used  to  compressively  reconstruct  the  signals  with  additive  Gaussian  noises  of 
[no  noise,  90,  30,  5]  SNRs.  The  LMMSE  optimizer  was  used  to  characterize  the  ID  imaging  performance  with  four 
different  operators  in  Fig.  7  and  the  LMMSE  optimization  is  compared  to  the  transpose  solution  of  the  operators  in 
Fig.  8.  The  LMMSE  optimization  is  generally  superior  to  the  transpose  solution  in  compressive  reconstruction. 
Looking  at  the  LMMSE  results,  the  ideal  operators  show  very  similar  normalized  root-mean-square  deviations 
(NRMSD)  with  the  VH  operators.  As  the  SNR  increases,  the 
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NRMSD  steeply  improves  in  general.  In  the  case  of  SNR=5,  the  NRMSDs  worsen  at  100.  In  the  numerical 
experiments,  four  system  operators  of  Hadamard,  PCA,  Random  Bernoulli,  and  Random  Gaussian  are  tested  with 
four  optimizers  of  LMMSE,  TV  minimization,  and  LI  minimization.  The  number  of  resolution  elements  are  64  for 
both  input  and  output  field-of  views. 

The  VH  implementation  is  not  much  inferior  to  the  ideal  system  operator  in  condition  number  and  log(MSE).  Now 
the  index  variation  Dn  is  small  enough  not  to  cause  a  significant  degradation.  The  results  are  promising,  showing 
that  the  performance  of  the  noise-free  VH  system  can  be  very  close  to  the  ideal  system.  The  x-axis  indicates  the 
sensing  rate  for  the  compressive  measurement.  So  100%  sensing  rate  means  64  measurements  for  64  object 
samples.  Since  optical  systems  cannot  measure  the  negative  values,  the  dual  rail  signal  was  used  subtracting  one 
positive  mask  with  the  other  positive  mask.  The  maximum  absolute  column  sum  normalization  is  also  taken  into 
account  to  preserve  the  photon  number  in  the  optical  measurement. 
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Fig.  7:  LMMSE  optimization  for  ID  imaging  operators:  (a)  Hadamard,  (b)  PCA,  (c)  Random 
code  with  Bernoulli,  and  (d)  Random  code  with  Gaussian  and  orthogonal  ization. 
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Fig.  8:  Transpose  for  ID  imaging  operators:  (a)  Hadamard,  (b)  PCA,  (c)  Random  code  with 
Bernoulli,  and  (d)  Random  code  with  Gaussian  and  orthogonalization. 


B.3  -  Task  3:  Investigate  the  use  of  completely  general  3D  index  distributions  for  the 
realization  of  non-identity  optical  mappings. 

This  task  has  been  pursued  in  two  thrusts.  The  first  thrust  has  leveraged  the  multi-linear  formulation  developed  in 
Task  1  to  realize  non-identity  optical  transforms.  Because  the  computational  complexity  of  the  resulting 
optimization  problem  prohibits  extension  to  high  dimensions  we  have  also  pursued  meta-material  solutions  to  this 
class  of  problem.  Once  again  in  accordance  with  kick-off  meeting  DARPA  feedback,  we  have  limited  our  attention 
to  identity  operators  in  this  second  thrust. 

B.3. a  Multi-Linear  Hadamard  Imager 

The  input/output  pairs  are  8x8  matrices.  Each  of  the  64  inputs  is  a  different  "point  source."  The  corresponding 
desired  output  matrix  is  the  output  of  the  operator 


-+1 

+1 

+1 

+1 

+1 

+1 

+1 

+1- 

+1 

+1 

+1 

+1 

-1 

-1 

-1 

-1 

+1 

+1 

-1 

-1 

-1 

-1 

+1 

+1 

+1 

+1 

-1 

-1 

+1 

+1 

-1 

-1 

+1 

-1 

-1 

+1 

+1 

-1 

-1 

+1 

+1 

-1 

-1 

+1 

-1 

+1 

+1 

-1 

+1 

-1 

+1 

-1 

-1 

+1 

-1 

+1 

-+1 

-1 

+1 

-1 

+1 

-1 

+1 

-1- 

16 

Approved  for  public  release:  distribution  is  unlimited. 


applied  to  the  point  source  input.  In  all  tests,  there  are  5  masks.  The  spacing  between  masks  is  1mm,  for  a  total 
propagation  distance  of  5mm.  The  oversampling  factor  is  28.  The  8x8  mask  pixels  are  spaced  10 fim  apart  in  both 
the  x  and  y  dimensions.  The  wavelength  is  .  5  \im .  1000  iterations  were  run,  and  two  tests  have  finished:  one  in 
which  only  the  center  mask  is  allowed  to  vary,  and  the  other  in  which  all  5  masks  are  allowed  to  vary. 

Figure  la  shows  the  intended  52nd  output;  Figure  lb  shows  the  51st  output  realized  in  the  test  with  the  single 
active  mask;  and  Figure  lc  shows  the  51st  output  realized  in  the  test  with  the  5  active  masks.  While  the  case  with 
the  5  active  masks  yields  an  output  closer  to  the  intended  one,  neither  test  realized  the  output  closely.  The 
average  RMSE  between  the  desired  output  and  the  realized,  normalized  output  is  plotted  in  figure  2. 


RMSE  v.  #  Masks 


#  Masks 

Figure  2:  Output  RMSE  versus  number  of  masks  for  the  Fladamard  index  distribution. 

B.3.b  Metamaterials  Re-Imager 

Quasi-conformal  transformation  optics,  or  QCTO,  was  invented  by  Li  and  Pendry  [1]  in  order  to  yield  designs  that 
were  more  manufacturable.  For  2-D  optical  systems,  in  which  light  only  propagates  in  two  dimensions,  QCTO  can 
be  used  to  produce  non-magnetic,  dielectric-only  material  prescriptions.  For  cylindrically-symmetric  optical 
systems,  QCTO  can  be  used  to  find  material  prescriptions  such  that  the  permeability  tensor  is  non-unity  only  in  the 
azimuth  direction.  These  prescriptions  are  more  likely  to  be  manufacturable  than  general  TO  prescriptions,  as 
most  magnetically-coupled  metamaterials  provide  a  magnetic  response  in  only  one  direction  [2]. 
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Material  prescriptions  that  have  a  narrower  range  of  constitutive  parameters  are  also  more  likely  to  be 
manufacturable.  This  paper  analyzes  the  effect  of  optimizing  the  shape  of  a  concave-piano  lens  designed  with 
QCTO  in  order  to  make  the  range  of  constitutive  parameters  as  narrow  as  possible.  QCTO  is  applied  to 
Maxwell's  fish-eye  lens,  reshaping  it  into  a  concave-piano  lens  that  can  image  a  curved  surface  onto  a  flat 
detector,  as  illustrated  in  Figure  1.  Various  parameters  of  the  shape  of  the  concave-piano  lens  are  varied,  and 
the  effect  on  the  range  of  constitutive  parameters  is  analyzed. 


Monocentric  Lens 


Figure  9:  Monocentric  Lens  Coupled  to  a  Concave-Piano  Reimager.  The  spherical  image  plane  of  the  monocentric  lens  is 

imaged  onto  a  flat  detector  plane. 


In  Figure  1,  the  concave-piano  lens  is  shown  reimaging  a  monocentric  lens.  Monocentric  lenses  are  spherically 
symmetric  lenses  that  have  a  hemispherical  image  surface.  Over  the  last  several  years,  several  promising  ways  to 
reimage  the  monocentric  lens'  curved  image  have  been  proposed  and  prototyped,  including  capturing  the  surface 
with  micro-camera  relays  [3]  and  transmitting  the  image  to  a  flat  detector  plane  via  fiber  bundles  [4].  The  concave- 
piano  reimager  here  would  represent  yet  another  method  of  reimaging  the  monocentric  lens. 


Figure  10:  The  relationship  between  the  curvilinear  quadrilaterals  and  the  transformation  .  This  figure  illustrates  the 
relationship  between  the  curvilinear  quadrilaterals  and  the  transformation. 


As  illustrated  in  Figure  2,  the  transformation  that  reshapes  Maxwell’s  fish-eye  lens,  a  sphere,  into  a  concave -piano 
shape  is  designated  T0^3dp\  U^3^  ->  U^3d\  U^3d^  and  U^3d^  are  three-dimensional  spaces  described  by  the  Cartesian 
coordinates  (x0,y0,z0)  and  (xcp,ycp,zcp)  respectively,  as  well  as  by  the  cylindrical  coordinates  (po,0o,zo)  and 
( Pep,  0Cp'zcp)  respectively.  The  z0-  and  zcp-  axes  are  the  two  optics’  axes  of  symmetry.  This  3-D  transformation  is 
related  to  the  2-D  conformal  transformation  T0^cp  \  U0  ->  Ucp  as  follows: 


( 


xcp>  fcp>  zcp 


) 


To(™i(xo,y0’Zo) 


Ocp’Pcp)  T0 


^cp 


o-^cp  Oo,  Po) 

4>o 


(30) 
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The  components  of  the  constitutive  parameter  tensors  then  are  specified  by  (  2  )-(  5  )  [2] : 

8 Pep  (ZCp>Pcp)  —  8Zcp  (Zcp>  Pep)  ~  £  (j'cp^ofecpi  Pcp)i^  * 


8  (^cp^ofecpf  Pcp)^) 

t(t>cpy.z,cp>Pcpj  -  r  tz  7  T\T 

| Jo^cp  y* cp^>o\zcp>  Pep) )\ 
Ppcp(Zcp>  Pep)  =  Pzcp(zcp>  Pcp)  =  1 


8<Pcp{.Zcp>  Pep) 


P<t>cp  (zcp>  Pep) 


f  Pep 


|/o->Cp  (j'cp^O  (zcp>  Pep))  |  P° 


y 


(31) 

(32) 

(33) 

(34) 


In  these  equations,  £pcp(zcp,  pcp),  £^cv{zcpt  pcp),  and  £Zcp(zcp,pcp)  are  the  components  of  the  permittivity  tensor  in 
the  new  optic,  oriented  in  the  radial,  azimuthal,  and  longitudinal  directions  respectively; 

Ppcp(zcp>Pcp)>  P(pcv  (zcp>Pcp)i  and  Pzcp  (zcp>  Pep)  are  the  corresponding  components  of  the  new  optic’s  permeability 
tensor;  e(ztp)  is  the  isotropic  permittivity  of  the  original  optic  (Maxwell’s  fish-eye);  p0  is  the  radial  coordinate  of 
the  preimage  of  ( zcp,pcp )  under  T0^cp  (z0,p0);  Tcp^0  =  T0^cp_1;  and  \j0^Cp(zo>  Po)\  denotes  the  determinant  of  the 
Jacobian  of  T0^cp(z0> p0 ).  Therefore,  in  the  material  described  by  (  2  )-(  5  ),  all  light  travelling  in  meridional  planes, 
regardless  of  polarization,  experiences  a  non-magnetic  response  with  index  of  refraction  described  by  (  6  ): 


ncp{_zcp>  Pep)  ~  ' 


n0  (^cp-»o  ( zcp>  Pcp)) 


J  \j  o—>cp  (z  cp>  Pcp)) | 

where  7i0  (z0,p0 )  is  the  original  index  distribution  in  Maxwell’s  fish-eye: 


n0  ( zo>  Po ) 


n0 


,  |  V  +Po 

Rn2 


(35) 


(36) 


where  n0  is  an  arbitrarily-chosen  scale  factor,  and  R0  is  the  radius  of  Maxwell’s  fish-eye,  which  may  be  arbitrarily 
set  to  1  mm. 


In  order  to  obtain  numerical  samples  of  the  constitutive  parameters  of  the  new,  reshaped  lens,  a  discrete  set  of  input 
points  and  corresponding  output  points  of  the  conformal  function  ( zcp,pcp )  =  T0^cp(z0,  p0)  are  first  found 
numerically,  by  way  of  the  following  methodology.  Initially,  two  “curvilinear  quadrilaterals,”  or  closed  chains  of 
four  arcs,  are  specified:  Q0  and  Qcp.  These  two  quadrilaterals  are  shown  in  Figure  2.  Q0  is  specified  in  terms  of  the 
radius  of  Maxwell’s  fish-eye,  R0 ,  which,  arbitrarily,  may  be  set  to  1  mm.  Q0  is  the  image  under  ^0(x0,y0,z0)  of  a 
closed  surface  that  mostly  aligns  with  the  outer  shell  of  Maxwell’s  fish-eye.  Here,  TTo  \  U^3^  ->  Ucp  designates  the 
projection  from  U^3d^  onto  UQ.  Similarly,  TTcp :  U^3d^  ->  Ucp  designates  the  projection  from  U^3^  onto  Ucp: 

TTo(x0,y0,z0)  =  (z0,Jx02  +y02)  =  Oo.Po)  (37) 

TpcpiXcp.ycp'Zcp)  \^cpi 
Hence  Qcp  sets  the  shape  of  the  reimager:  it  is  the  image  of  the  boundary  surface  of  the  reimager  under 
J?cp(xcp>ycp>zcp)-  The  left  arc  of  Qcp  corresponds  to  the  concave  front  surface.  The  right  line  segment  corresponds 
to  the  flat  back  detector  surface.  Then,  with  numerical  grid-generation  techniques,  as  described  in  [5]  and  [6],  grid- 
point  arrays,  (z0[i,j],p0[i,j]  )  and  (zcp[i,j],pcp[i,j]  ),  are  computed  on  and  inside  Q0  and  Qcp  respectively,  such 
that: 


" cp 


T  Tcp 


)  (zcp>  Pcp) 


(38) 


^o^cp  (zo  U>j]>  Po  [Uj])  =  (zcp[i>j],Pcp[i>j])  (  39) 

where  J^cp(z0,p0)  is  conformal. 

Then,  for  each  point  (zcp[i,j],pcp[i,j]),  an  approximation  of  ncp(zcp[i,  j],pcp\i,j])  is  computed  with  the  equation: 
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X 


(40) 


where: 


and 


*  j,  = _ no _  I  l^° I u 

ncP  i’j  t  +  (z0[i,j]/R0y  +  ( p0[i,j]/R0y  x  |;^cP|iy 


|/ P->0  |  j  J 


vz0[i  +  1  ,j]~z0[i,j] 

z0[i  +  1J]  -z0[i.;]l 

A  T] 

PoU  +  1  ,j]~Po[i.j] 

PoU  +  1  J\-Po[i.j\ 

i  U( 

A  T]  \ 

1/ 


^cp  I  i 


[zcp[i  +  l,j]  -  zcp[i,j] 

zcp[i  +  1J]  -ZcP[Uj]-\ 

A  T] 

PcP[i  +  l,j]~  pcp[i,j] 

PcpU  +  1J]  -  PcP[iJ ] 

Ap  \ 

are  the  first  order  finite  difference  numeric  approximations  of  |/^0(^,^(Zo[b /LPo  [b7]))|  and 
(^cp^f^cp  [i.jl  Pco [i.;]))  I  respectively. 


(41) 


(42) 


As  long  as  the  quadrilaterals’  lateral  sides  are  aligned  with  the  optics’  object  and  image  surfaces,  arbitrary 
quadrilaterals  accomplish  the  design  goal.  Different  quadrilaterals  Q0  and  Qcp,  however,  yield  different  constitutive 
parameters.  The  effect  of  several  different  parameters  of  the  quadrilaterals  on  the  range  of  the  refractive  index  is 
thus  assessed  with  the  following  procedure.  First,  a  particular  Q0  is  chosen.  Then  Qcp  is  defined  by  a  parameterized 
formula.  The  value  of  each  Qcp  parameter  is  varied  over  several  test  trials,  while  the  others  are  held  fixed.  For  each 
test  case,  the  metric  An™p  is  calculated,  as  defined  in  (  14  ): 


max  { ncp[i,j ]} 


A  nw  = 

CP  {ncp[i>j]} 


(43) 


Dn  in  (  14  )  denotes  the  set  ofpoints  (zcp[t,y],pcp [(,)])  £  fl^60^  c  Ucp  is  the  set  ofpoints  on  and 

interior  to  Qcp  to  which  Jj>cp  (xcp,ycp,zcp),  defined  in  (  9  ),  projects  the  points  inside  the  reimager  needed  to 
reimage  a  ±60°  field.  A  fi™p  is  therefore  the  ratio  of  the  maximum  index  to  the  minimum  index.  It  is  useful  because  it 
does  not  depend  on  the  arbitrarily-chosen  scale  factor  n0  in  (  7  )  nor  the  arbitrarily-chosen  radius  of  Maxwell’s  fish- 
eye  R0.  Finally,  the  relationship  between  the  Qcp  parameters  and  A n™p  is  used  to  assess  the  effects  of  the 
quadrilateral  parameters  on  index  range  and  determine  the  final  design. 

The  chosen  structure  of  Q0  is  diagrammed  in  Figure  3.  The  bottom  side  of  Q0  is  the  diameter  of  Maxwell’s  fish-eye 
lens  that  is  on  the  optical  axis:  {(z0, 0)  E  U0\  \z0  \  <  R0}.  Meanwhile,  the  left  and  right  sides  are: 

{( +R0  cos  6  ,  R0  sin  O')  E  U0  \  0  <  0  <  cp0 },  where  cp0  =  80°;  where  the  negative  sign  corresponds  to  the  left  side; 
and  where  the  positive  sign  to  the  right  side.  The  top  side  is  symmetric  in  z0.  Each  half  of  this  symmetric  top  side  is 
defined  by  the  cubic  parametric  equations  in  (  15  )  and  (16)for0<s<l,  where  the  top  signs  in  the  t=F’  and  ‘±’ 
correspond  to  the  left  part;  and  the  bottom  signs  to  the  right  segment.  (  15  )  specifies  zoT(s)/R0 ,  the  ratio  of  the  z0- 
coordinate  of  the  top  side  to  R0,  as  a  function  of  the  parameter  s.  (  16  )  similarly  specifies  poT(s)/R0 ,  the  ratio  of 
the  p0  -coordinate  to  R0  as  a  function  of  s.  The  top  is  configured  such  that  it  forms  a  90°  angle  with  the  side,  such 
that  it  is  level  at  the  midpoint,  at  height  h0  =  .8 R0,  and,  arbitrarily,  such  that  |zo/(0)  |  =  3  x  R0  cos  cp0. 

Z0T(s)/R0  =  +  COS  <Po±\  COS  <Po  5  +  j±  ^  COS  <p0  —  ij  s2  +  [+1.75  cos  <p0  +i]s3 

Paris) /R0  =  sin  cp0  —  isin  (p0  s  + 1 —  —  sin  (p0  +  3/i0js2  +  {1.75  sin  cp0  —  2h0}s3 


(44) 

(45) 
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While  it  might  be  useful  in  general  to  study  variations  in  Q0,  in  addition  to  Qcp,  this  analysis  only  varies  Qcp 
parameters.  For  all  tests,  <p0  =  80°  and  h0  =  .8 R0. 

For  this  analysis,  Qcp  is  defined  with  the  parameters  Rcp ,  W0 ,  0,  HD,  and  a  depicted  in  Figure  4.  Rcp  is  the  radius  of 
the  monocentric  imager.  Since  this  reimager  is  designed  based  on  the  monocentric  lens  of  [7],  Rcp  =  11.9885  mm. 
W0,  meanwhile,  is  the  length  of  the  reimager  along  the  optical  axis.  As  shown  in  Figure  4,  the  bottom  side  of  Qcp  is 
a  horizontal  line  segment  on  the  zcp-axis:  {(zcp,  0)  E  Ucp\Rcp  <  zcp  <  Rcp  +  W0).  The  left  side,  additionally,  is  a 
subset  of  the  image,  under  TTcp(xcp,ycp,zcp ),  of  the  monocentric’s  image  surface.  It  subtends  angle  The 
monocentric  lens’  center  of  curvature  corresponds  to  the  origin  of  Ucp  \  hence  the  left  side  is: 

|(ftcp  cos  6  ,  Rcp  sin  9 )  E  Ucp  \  0  <  6  <  jJ.  ^  must  be  at  least  the  desired  half  field  of  view.  By  contrast,  the  right 
side  is  a  straight  vertical  line  segment  of  length  where  HD  is  the  hill  diameter  of  the  flat  back  surface,  where  the 
detector  can  be  placed.  It  is  the  set  |(ftcp  +  W0 ,pcp)  £  Ucp :  0  <  pcp  <  ^|. 


P  ^rn 

Kcp  cp 


Figure  12:Qcp.  Qcp  is  parameterized  by  Rcp,  W0,  0,  ^  and  a  (a  is  directly  related  to  H) 


The  top  side  of  Qcp ,  shown  in  red  in  Figure  4,  is  defined  by  parametric  equations,  (z£pop^(s),  p^pp\s)^,  which  are 

polynomial  functions  of  s  for  0  <  s  <  1.  For  smoothness,  the  curve  is  chosen  to  have  third-order  geometric  (G3) 
continuity  -  i.e.  if  5  represented  time,  the  path’s  linear  and  angular  acceleration  would  be  continuous.  According  to 
[8],  the  set  of  polynomials  of  the  lowest  order  that  can  be  parametric  equations  of  a  G3  curve,  given  any  particular 
set  of  starting  and  ending  points;  starting  and  ending  tangential  angles;  starting  and  ending  curvatures;  and  starting 
and  ending  angular  accelerations,  are  the  seventh-order  ones  that  [8]  describes,  in  which  the  coefficients  are 
functions  of  these  eight  boundary  conditions  along  with  six  other  parameters:  9^  Rs?  and  ^6-  The  top  side 
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of  Qcp  is  thus  defined  by  such  a  pair  of  polynomials.  The  starting  point  (z£p ,  )  and  ending  point  {z^p  ,  p^)  used 

are  the  uppermost  points  of  the  left  and  right  sides  respectively: 

(4p'  Pep  )  =  (Rcp  cos  | ,  Rcp  sin  y)  f  46  ) 

(?cp\pcj)=  (Rcp  +  Wo.y)  (47) 

Because  the  boundaries  of  Qcp  are  isolines  of  the  conformal  Tcp^(zcp,  pcpf  moreover,  the  right  and  left  sides  of 
Qcp  necessarily  meet  the  top  and  bottom  sides  of  Qcp  at  right  angles.  Therefore,  the  starting  and  ending  tangential 
angles  of  the  top  curve  are  the  ones  that  form  right  angles  with  the  adjacent  sides.  The  starting  tangential  angle  is 
\pi=^  and  the  ending  one  is  ipf  =  0°.  In  all  test  cases,  the  initial  and  final  angular  “acceleration”  (second-order 
derivative  with  respect  to  5)  is  arbitrarily  set  to  0:  kt  =  0,  kf  =  0.  Additionally,  the  initial  curvature  is  set  to - , 

while  the  final  curvature  is  set  to  0:  Kt  =  —  1  /Rcp,  Kf  =  0.  ”G3-Splines  for  the  Path  Planning  of  Wheeled  Mobile 
Robots”  (  [8])  furthermore  states  that,  given  the  boundary  conditions,  an  approximately  optimally  smooth  pair  of 
polynomials  is,  in  many  “typical”  cases,  the  one  in  which  p1  and  p2  are  both  set  to  the  distance  between  the  starting 
and  ending  points,  and  the  other  p’s  are  set  to  0: 


Vl=ri2  =  II  (Zcp.Pcp)  -  (Zcp’Pcp)  II 

(48) 

C3 

00 

II 

II 

on 

II 

Gs 

II 

O 

(49) 

Based  on  experimentation,  we  find  that  if 

Vl=V2=  4 \(zcp-Pcp)  -  (Zcp’Pcp)W 

(50) 

then  variation  of  the  parameter  a  linearly  affects  the  overall  height  H  of  the  reimager.  is  illustrated  in  Figure  4.) 

Thus,  ZcP°p\s )  and  p^pv\s')  are  provided  by  (  22  )  and  (  23  ).  The  a  and  b  coefficients  in  these  formulas  are 

obtained  by  fitting. 

zcP°p)(s)  =  a0  +  a±s  +  a2s2  +  a3s 3  +  a4s 4  +  asss  +  a6s6  +  a7s7  (51 ) 

Pcp0P\s )  =  b0  +  bxs  +  b2s 2  +  b3s3  +  b4s 4  +  bsss  +  b6s6  +  b7s7  ( 52  ) 


Table  1  lists  the  tests  performed.  The  first  test  is  intended  to  determine  the  effect  on  index  range  of  W0,  the  depth  of 
the  reimager.  Test  2  varies  HD,  in  order  to  assess  the  effect  of  the  height  of  the  back  detector  plane  on  the  index 
range.  Finally,  the  third  test  varies  the  parameter  a  defined  in  (  21  ).  This  test  effectively  varies  the  overall  height  H 
of  the  reimager  with  each  trial,  and  thus  assesses  the  relationship  between  H  and  index  range. 

Table  3:  List  of  Tests 


Test  # 

Qcp  Parameter 

W0/Rcp 

cp 

Hd/RcV 

a 

1 

{1.75,2,  ...,5.5} 

180° 

2.5 

1 

2 

3.5 

180° 

{.75, 1, ...  ,5.5} 

1 

3 

2.5 

180° 

2.5 

{0.5, 0.6,  ...,1.5} 

Width  test 

Using  the  parameters  listed  in  the  first  row  of  Table  1,  the  first  test  investigates  the  effect  of  W0  on  index  range.  W0 
is  varied  from  1.7SRcp  to  5.5 Rcp  in  increments  of .  25 Rcp.  The  resulting  A fVfp  measurements  are  plotted  in  Figure  5. 
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15 


An^  v.  WQ 


A n™v  decreases  with  W0.  The  reason  for  this  decrease  is  found  by  analyzing  Figure  6a  and  b.  Figure  6a  and  b 
illustrate  the  paths  of  ray  bundles  through  the  reimager  for  the  two  cases  in  which  W0  =  1.7SRcp  and  W0  =  5.5 Rcp 
respectively.  The  ray  bundles  in  the  longer  reimager  are  bent  more  gradually,  which  necessitates  a  narrower  index 
range,  in  agreement  with  the  finding  that  a  longer  reimager  corresponds  to  a  lower  A n™p.  Based  on  these  results,  it  is 
concluded  that  a  longer  reimager  is  more  likely  to  be  manufacturable.  The  larger  W0,  the  narrower  the  range  of 
index  values  is  required. 


=  1.75  R 

Q  cp 


-2  - ■ - ■ - ■ 

■0  1  2  3 


z.  R 

cp  cp 

Figure  14a:  Ray  Bundles  through  Reimager  with  W0  = 
1.75  Rcp 


Wn  =  5.5  R 

0  cp 


0  2  4  6 


Figure  6b:  Ray  Bundles  through  Reimager  with  W0  =  5.5 Rcp 


Back  height  test 

In  the  second  test,  the  effect  of  the  height  of  the  back  surface  of  the  reimager,  or  image  height,  on  the  index  range  is 
gauged.  Qcp  is  generated  with  the  parameters  in  the  second  row  of  Table  1.  HD,  the  full  height  of  the  back  plane  of 
the  reimager,  is  varied  from  .  75 Rcp  to  5.5 Rcp  in  increments  of .  25 Rcp.  The  different  Qcp  s  corresponding  to  each 
test  case  are  shown  in  Figure  7.  The  test  is  created  with  the  hypothesis  that  the  index  range  is  likely  to  be  narrower 
when  Hd  is  close  to  values  that  correspond  to  an  average  unit  magnification  of  the  image:  HD  «  nRcp. 
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Qcp  for  Different  HD 


Figure  15:  How  the  Top  Side  of  Qcp  varies  with  HD.  The  top  side  of  Qcp  terminates  at  the  coordinate  pcp  =  As  HD  is 
increased,  the  top  side  of  Qcp  becomes  smoother.  It  also  becomes  taller,  which  increases  the  conformal  module.  As  HD  ->  0, 
the  backside  of  Qcp  becomes  increasingly  compressed.  The  case  where  HD  =  2 Rcp  is  emboldened,  since  it  has  the 

narrowest  index  range. 

The  computed  A ncp  for  each  value  of  HD  is  plotted  in  Figure  8.  For  HD  <  2 Rcp,  A iVfp  decreases  with  Rcp.  The 
narrowest  index  range  occurs  for  HD  =  2 Rcp.  Then,  as  HD  is  increased  beyond  2 Rcp,  the  index  range  widens. 


1  2  3  4  5 


^ D  /^cp 

Figure  16:  A fVfp  v.  HD.  A n™p  decreases  with  HD  for  HD  <  2 Rcp  and  increases  with  HD  for  HD  >  2 Rcp.  (The points 
with  the  red  x  ’s  represent  cases  in  which  Qcp  does  not  encompass  the  area  needed  for  all  the  rays  paths 

originating  within  60°  of  the  axis  of  symmetry.) 


This  trend  is  explained  by  Figure  9a,  b,  and  c,  which  show  traces  of  ray  bundles,  from  field  angles  of  0°  and  60°, 
through  the  reimager  for  the  cases  in  which  HD  =  0.75 Rcp,  HD  =  2 Rcp,  and  HD  =  5.5 Rcp  respectively.  (The  60° 
bundle  in  the  HD  =  0.75 Rcp  case  is  not  shown  because  this  Qcp  does  not  encompass  enough  points  to  image  a  60° 
field  of  view.)  On  the  one  hand,  in  the  case  where  HD  =  0.7SRcp,  the  rays  must  be  abruptly  straightened  in  order  to 
fit  inside  the  compressed  back  of  the  reimager.  This  abrupt  change  in  direction  necessitates  a  wider  index  range.  On 
the  other  hand,  in  the  case  where  HD  =  S.SRcp,  the  60°  ray  bundle  is  elongated,  compared  to  the  one  in  the  case 
where  HD  =  2 Rcp.  A  wider  index  range  is  also  needed  for  this  stretching.  In  this  way,  these  ray  bundle  diagrams  in 
Figure  9a-c  provide  a  physical  interpretation  of  the  trends  in  Figure  8. 
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Figure  1 7a:  Ray  Bundles  for  the  Case  in  FiSure  9b:  Rav  Bundles  for  the  Case  Figure  9c:  Ray  Bundles  for  the  Case  in 


Which  Hd  =  0.75 Rcp  in  Which  H D  =  2 Rcp  Which  H D  =  5.5 Rcp 

Traces  of  ray  bundles,  from  field  angles  of  0°  and  60°.  (The  60°  bundle  in  the  HD  =  0.75  Rcp  case  is  not  shown  because  this 
Qcp  does  not  encompass  enough  points  to  image  a  60°  field  of  view.)  In  the  case  where  HD  =  0.7SRcp,  the  rays  must  be 
abruptly  straightened  in  order  to  fit  inside  the  compressed  back  of  the  reimager.  In  the  case  where  HD  =  5.5  Rcp,  the  60°  ray 
bundle  is  elongated,  compared  to  the  one  in  the  case  where  HD  =  2 Rcp.  Thus  the  HD  =  0.75 Rcp  and  HD  =  5.5 Rcp  cases 
correspond  to  a  higher  A ncp  and  Anfp  than  the  case  in  which  HD  =  2 Rcp. 


This  test  suggests  that  there  is  an  optimal  HD  for  the  reimager.  In  this  test,  the  optimal  value  is  HD  =  2 Rcp.  It  is 
anticipated  in  general,  however,  that  the  optimum  might  depend  on  the  other  Qcp  parameters,  since  the  value  2  Rcp 
corresponds  to  the  point  at  which  the  set  of  Qcp  repositions  the  minimum  \j^CpCTcp^(zcp> Pcp))\- 


Overall  height  test 

The  third  test  seeks  to  determine  the  effect  of  the  overall  height  of  the  reimager  on  index  range.  It  uses  the  Qcp 
parameters  listed  in  the  third  row  of  Table  1.  In  order  to  make  the  overall  height  H  vary,  the  parameter  a  in  (  21  )  is 
changed  from  .  5  to  1.5  in  increments  of .  1.  The  overall  height  of  Qcp  increases  linearly  with  a.  This  relationship  is 
demonstrated  in  Figure  10,  which  illustrates  the  Qcp  corresponding  to  each  tested  a,  as  well  as  by  the  measurements 
of  overall  height  H  as  a  function  of  a  plotted  in  Figure  1 1 .  The  test  is  undertaken  with  the  expectation  that  the  Qcp  ’s 
with  values  of  a  nearest  1  yield  the  narrowest  index  ranges,  since  their  top  sides  are  likely  to  be  the  smoothest, 
according  to  [8]. 


2.5 
2 

1.5 

Pcp/^cp  ^ 

0.5 

0 

-10  1  2  3  4 


Qcp  f°r  Different  a  Values 


\ 

a  - 1.3 

a  =  1 
a  =  0  5 


Figure  18:  Qcp  for  Different  Values  of  a.  The  parameter  a  is  used  to  vary  the  overall  height  of  the  reimager. 
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a 


Figure  19:  Hv.  a.  For  each  test  case,  the  overall  height  H  of  the  reimager  is  measured.  (The  parameter  H  is  illustrated  in 

Figure  4.)  H  varies  linearly  with  a. 

The  calculated  A n™v  for  each  test  are  plotted  in  Figure  12  respectively.  A fVfv  is  minimal  around  a  =  .7. 


Figure  20:  A fVfp  v.  a.  The  minimum  A fVfp  is  at  a  ~  .7. 

To  help  explain  the  trend  in  Figure  12,  Figure  13a,  b,  and  c  show  traces  of  ray  bundles,  from  field  angles  of  0°  and 
60°,  through  the  reimager  for  the  cases  in  which  H  =  2.9 Rcp,  H  =  3.2 Rcp,  and  H  =  4.17 Rcp  respectively.  The 
difference  between  the  ray  bundle  trajectories  in  the  H  =  2.9 Rcp  case  and  the  H  =  3.2 Rcp  case  does  not  appear 
significant.  On  the  one  hand,  in  the  case  where  H  =  4.17 Rcp,  the  rays  must  arch  higher.  This  longer,  more  arched 
path  explains  the  necessity  for  the  wider  index  range  associated  with  the  greater  heights.  In  this  way,  these  ray 
bundle  diagrams  in  Figure  13a-c  provide  a  physical  interpretation  of  the  trends  in  Figure  12. 
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Figure  21a:  Ray  Bundles  for  the  Case  Figure  13b:  Ray  Bundles  for  the  Case  Figure  13c:  Ray  Bundles  for  the  Case 

in  Which  H  =  2.9 Rcp  in  Which  H  =  3.2 Rcp  in  Which  H  =  4.17 Rcp 

Traces  of  ray  bundles,  from  field  angles  of  0°  and  60°,  through  the  reimager  for  the  cases  in  which  (a)  H  =  2.9  Rcp,  (b) 

H  =  3.2 Rcp,  and  (c)  H  =  4.17 Rcp.  In  the  case  where  H  =  4.1 7Rcp,  the  rays  must  arch  higher.  This  longer,  more  arched 

path  explains  the  necessity  for  the  wider  index  range. 

The  test  results  reveal  that  a  and  thus  H  have  an  optimal  value.  More  vertical  room  for  the  reimager  beyond  the 
height  associated  with  this  optimum  is  not  seen  to  improve  the  manufacturability  -  in  fact  requiring  the  light 
trajectories  to  take  paths  that  are  too  high  necessitates  greater  index  variation. 

After  studying  the  effect  of  different  reimager  shape  parameters  on  the  range  of  the  constitutive  parameters,  the 
material  prescriptions  for  a  reimager  that  can  be  coupled  to  the  particular  monocentric  lens  specified  by  [7]  are 
provided,  and  the  optical  system  is  ray-traced  to  verify  that  it  images  as  intended.  The  reimager  is  designed  using  a 
Qcp  with  W0  =  4 Rcp,  0  =  180°,  HD  =  3.5 Rcp,  and  H  =  4.549 Rcp.  The  calculated  constitutive  parameters  for  the 
reimager  needed  to  reimage  a  ±60°  field  are  were  computed.  Using  these  parameters  bundles  of  meridional  rays  are 
traced  through  a  cross-section  of  the  reimager.  Each  bundle  spans  ±25°,  approximately  the  angular  extent  of  the 
cones  of  light  entering  the  reimager  from  the  F/l  monocentric  lens.  (The  bundles  at  larger  field  angles  are  expected 
to  have  a  narrower  span,  due  to  surface  refraction.)  The  ray  trace  is  shown  in  Figure  15,  and  the  spot  sizes  are  shown 
in  Figure  16.  The  spot  sizes  of  the  30°  and  60°  field  angles  at  the  back  surface  are  less  than  2.5  microns  - 
comparable  to  the  size  of  the  fibers  bundles.  The  on-axis  rays  have  a  spot  size  of  less  than  6  microns.  According  to 
transformation  optics  theory,  since  the  points  imaged  by  Maxwell’s  fish-eye  lens  have  no  finite  extent,  at  least  in  the 
geometric  limit,  then  these  spot  sizes  should  also  be  0.  It  is  therefore  not  unreasonable  to  suppose  that  these  small 
measured  spot  sizes  might  be  a  precision  error,  either  in  the  discrete  model  of  the  constitutive  parameters,  or  in  the 
ray-tracer. 


60  = 
30  = 
0" 


z  ,  mm 

Figure  22:  Ray  Trace.  Bundles  of  meridional  rays,  from  field  angles  of0°,  30°  and  60°  are  traced  through  a  cross-section 

of  the  reimager.  Each  bundle  spans  ±25°. 
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Figure  23:  Spot  Sizes  of  Meridional  Rays 


B.4  -  Task  4:  Implications  for  MGRIN. 

The  MRGIN  performers  are  generally  unable  to  create  arbitrary  3D  refractive  index  distributions.  It  is  typical  for  a 
specific  manufacturing  process  to  yield  index  distributions  with  a  characteristic  symmetry  (e.g.,  layer  materials, 
rotationally  symmetric  materials,  etc.).  In  this  task  we  have  developed  a  practical  design/analysis  methodology 
that  can  capture  these  symmetries  and  enable  efficiency  optimization. 

B.4.a.  Refractive  index  basis  functions  as  a  compact  representation  for  complex  distributions. 

Recent  work  has  demonstrated  the  principle  of  transformation  optics  (TO)  [1-4]  as  a  method  to  design  gradient 
index  (GRIN)  optical  elements.  In  the  TO  design  method  the  desired  light  ray  trajectories  are  first  traced  through  a 
given  material.  A  transformation  is  then  applied  that  yields  the  necessary  spatial  distribution  in  permittivity  and 
permeability  of  that  material  to  yield  the  desired  ray  trajectories.  Unfortunately,  the  result  of  this  approach 
typically  yields  materials  which  are  highly  anisotropic,  and/or  magnetic,  and  thus  are  not  readily  fabricated.  To 
overcome  this  barrier  a  quasi-conformal  transform  (QCT)  optimization  approach  was  developed  that  removes  the 
need  for  magnetic  materials  and  minimizes  the  anisotropy  of  the  TO  design  [5].  The  resulting  optical  elements  only 
feature  spatial  variations  in  the  permittivity  (refractive  index).  While  the  QCT  approach  is  highly  effective  in  some 
cases  it  does  not  always  achieve  the  same  performance  as  the  original  TO  design  with  magnetic  and  anisotropic 
materials.  Thus  there  is  a  need  for  further  optimization  to  recover  the  performance  losses  that  come  from  TO/QCT 
GRIN  lens  designs. 

Optimization  techniques  for  GRIN  lenses  available  in  commercial  ray  tracing  software  (e.g.  Zemax,  CODE  V) 
typically  assume  a  functional  form  for  the  3D  spatial  index  distribution.  Thus  the  GRIN  lens  can  be  parameterized 
by  a  relatively  few  number  of  coefficients  and  easily  optimized.  The  GRIN  lens  designs  resulting  form  TO/QCT 
generally  cannot  be  represented  by  the  analytic  functional  forms  available  in  commercial  ray  tracing  software,  and 
thus  cannot  be  easily  optimized  using  these  packages.  Even  if  the  index  profile  is  manually  incorporated  into  the 
software  using  a  macro  or  add-on,  the  number  of  free  variables  (i.e.  voxels  within  the  GRIN  volume)  is  so  high  as  to 
make  optimization  infeasible.  The  optimization  issues  are  further  exacerbated  by  the  fact  that  algorithms  in  ray 
tracing  packages  often  heavily  rely  on  symmetry  and  assumptions  that  may  not  be  valid  with  these  GRIN  elements. 
Thus  there  is  a  need  for  an  efficient  way  to  represent  and  optimize  TO/QCT  GRIN  lens  designs. 

To  accomplish  this  goal,  we  propose  the  use  of  basis  functions  to  compactly  represent  the  3D  index  distribution  of 
the  GRIN  lens  by  a  relatively  few  number  of  coefficients.  Using  this  approach,  the  10A6-10A8  voxels  for  a  millimeter 
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scale  GRIN  lens  can  be  accurately  represented  by  less  than  1000  coefficients.  These  coefficients  can  then  be  used 
as  the  free  variables  in  an  optimization  algorithm  to  retrieve  the  performance  of  the  original  TO  design.  While  not 
as  compact  as  closed  form  purely  analytic  GRIN  lenses,  the  basis  function  approach  still  yields  a  feasible 
optimization  problem  and  gives  added  freedom  to  search  the  design  space. 

Some  relevant  questions  surrounding  this  approach  include: 

•  What  is  the  optimal  basis  for  defining  GRIN  optics? 

•  How  does  one  choose  which  and  how  many  coefficients  to  include  in  the  optimization? 

•  What  is  the  minimum  number  of  ray  and/or  field  points  necessary  to  properly  sample  the  space  within 
the  GRIN  lens? 

The  goal  of  this  initial  investigation  was  to  answer  some  these  questions  and  develop  heuristics  for  the  design  and 
optimization  of  the  GRIN  lenses  using  the  basis  function  approach. 

B.4.b.  A  novel  ray-tracing  engine  that  can  operate  on  these  index  representations 

A  novel  design  engine  was  developed  in  an  attempt  to  answer  these  questions.  The  GRIN  ray  tracing  algorithm  was 
based  on  the  work  for  Sharma  et  al.  [6].  This  ray  tracer  was  used  as  the  forward  model  in  the  design  engine,  and 
the  Covariance  Matrix  Adaptation  -  Evolutionary  Strategy  (CMA-ES)  algorithm  was  used  as  the  optimizer  [7].  The 
algorithmic  flow  through  the  design  engine  can  be  described  as  follows: 

1.  Input  an  initial  GRIN  lens  design  that  is  represented  within  a  given  basis 

2.  A  finite  number  of  basis  coefficients  are  kept  in  order  to  ensure  an  accurate  representation,  and  the  rest 
are  discarded  (i.e.  set  to  zero) 

3.  The  GRIN  lens  is  reconstructed  from  these  finite  coefficients,  and  several  bundles  of  rays  are  traced 
through  the  structure  and  a  figure  of  merit  is  calculated  (typically  RMS  spot  size) 

4.  The  CMA-ES  algorithm  adjusts  the  coefficients  and  the  rays  are  re-traced  through  the  updated  GRIN  lens 

5.  The  figure  of  merit  is  recalculated  and  kept  if  it  is  smaller  than  the  current  minimum 

6.  Steps  4  and  5  are  repeated  until  the  desired  FOM  is  achieved  or  a  fixed  number  of  iterations  has  been 
reached 

The  first  study  used  a  standard  GRIN  lens  with  a  quadratic  profile  in  the  transverse  dimension  and  a  constant  index 
in  the  axial  dimension  as  its  starting  design.  This  lens  was  chosen  as  it  is  well  known  in  the  literature  and  would 
serve  as  a  good  test  case  to  understand  the  workings  and  limitations  of  the  design  engine.  Additionally,  it  can  be 
easily  represented  in  Zemax  and  thus  provide  validation  for  our  forward  model.  The  input  rays  were  configured 
such  that  the  lens  was  focusing  collimated  rays  from  five  different  field  points  (0°,  ±15.35°,  ±30.7°)  at  the  back  face 
of  the  lens.  Two  different  basis  functions  were  used:  the  discrete  cosine  transform  (DCT)  and  Daubchies  wavelets. 

In  the  first  experiment,  only  the  largest  valued  coefficients  from  the  initial  fit  were  kept  for  optimization  as  they 
yielded  the  most  accurate  representation  of  the  initial  design.  The  wavelet  basis  optimized  effectively,  however 
the  DCT  basis  was  unable  to  improve  upon  the  initial  design  as  seen  in  Fig.  1.  The  reason  for  this  is  that  the  largest 
coefficients  of  the  initial  fit  were  only  along  the  transverse  dimension,  and  thus  the  design  engine  was  unable  to 
alter  the  axial  index  profile.  This  experiment  demonstrated  an  important  and  not  immediately  obvious 
consideration  when  attempting  to  optimize  GRIN  lenses  using  basis  functions:  The  coefficients  which  best 
represent  the  initial  design  do  are  not  necessarily  the  best  coefficients  for  optimization.  The  DCT  optimization  was 
repeated  using  low  frequency  coefficients  along  both  dimensions  and  proceeded  to  optimize  similarly  to  the 
wavelets,  even  though  the  initial  fit  was  somewhat  poor. 
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Figure  1 


The  second  experiment  attempted  to  answer  the  question  as  to  the  desired  number  of  coefficients  need  for 
optimization.  A  cursory  hypothesis  is  that  more  coefficients  would  give  the  optimizer  more  degrees  of  freedom 
and  lead  to  better  results,  but  at  the  cost  longer  computation  times.  To  test  this  hypothesis  optimization  for  the 
quadratic  GRIN  lens  were  repeated  with  an  increasing  number  of  coefficients  up  to  600.  It  was  found  that  when 
less  than  50  coefficients  were  used,  the  initial  reconstruction  of  the  GRIN  lens  was  so  poor  that  the  rays  were  not 
accurately  traced  through  the  lens  and  the  optimizer  was  unable  to  improve  on  the  design.  This  provides  initial 
insight  on  the  lower  bound  on  the  number  coefficients  needed  for  successful  optimization.  As  the  number  of 
coefficients  was  increased  the  FOM  reached  a  minimum  and  then  steadily  began  to  increase  as  shown  in  Fig.  2. 
This  was  true  for  both  basis  sets  with  the  optimal  number  of  coefficient  being  near  100.  Add  more  coefficients 
appears  to  have  created  too  large  of  a  solution  space  efficient  optimization.  The  wavelet  basis  achieved  a  slightly 
lower  figure  of  merit,  but  more  rapidly  diverged  from  the  optimal  result  when  too  many  coefficients  were  used. 


Quadratic  GRIN 


The  third  experiment  in  the  initial  study  attempted  to  determine  the  minimum  number  of  rays  needed  to 
accurately  determine  the  figure  of  merit.  Each  ray  traced  through  the  system  adds  to  the  optimization  time,  and 
thus  it  is  imperative  to  minimize  them.  However,  using  too  few  rays  leads  to  inaccurate  sampling  the  of  the  lens 
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and  poor  designs.  Thus  the  number  of  rays  needed  to  effectively  optimize  was  defined  as  the  minimum  number  of 
rays  such  that  if  the  ray  density  of  the  bundle  of  the  final  design  is  arbitrarily  increased  the  RMS  spot  size  does  not 
substantially  increase.  Table  1  shows  the  results  of  this  experiment  for  both  basis  functions.  It  was  found  that  using 
25  rays  per  bundle  in  the  DCT  basis  and  30  rays  per  bundle  in  the  wavelet  basis  were  required  to  produce  less  than 
a  1  pm  change  in  the  RMS  spot  size  when  300  rays  were  traced  in  the  same  bundle.  This  is  a  logical  result  as  the 
wavelet  basis  results  in  more  local  high  frequency  oscillations  than  the  DCT  which  needs  denser  sampling.  So  even 
though  the  initial  analysis  suggests  that  the  wavelet  basis  yielded  better  results,  it  required  more  rays  per  bundle 
(and  thus  a  slower  optimization  time)  to  achieve  that  result. 


Average  Increase  in  RMS  Spot  Size 
for  300  ravs  traced  per  bundle  (um) 


Number  of  Rays 
per  Bundle  Used 
for  Optimization 

DCT  Fit 

Wavelet  Fit 

7 

4.4 

54.2 

15 

2.8 

7.7 

20 

1.2 

3.3 

25 

0.7 

1.6 

30 

0.5 

0.9 

Table  1 


Combining  the  results  described  thus  far  one  can  generally  conclude  that  the  basis  function  approach  to  GRIN  lens 
design  needs  a  sufficient  number  of  coefficients  to  effectively  represent  the  initial  design,  but  not  so  many 
coefficients  as  to  create  an  unsearchable  space.  Furthermore,  the  distribution  of  those  coefficients  should  be 
chosen  to  ensure  that  space  is  not  over  constrained  by  the  geometry  starting  design.  Finally,  the  number  of  rays 
per  bundle  should  be  dense  enough  as  to  effectively  sample  the  index  distribution  within  the  lens.  Applying  these 
heuristics  to  the  quadratic  GRIN  singlet  the  basis  function  approach  was  able  to  improve  RMS  spot  size  over  the 
full  field  by  an  order  of  magnitude  and  approach  diffraction  limited  resolution,  as  shown  in  Table  2. 


Field 

Point 

Zemax  Optimized 
RMS  Spot  Size  (pm) 

Optimized  DCT  Fit 
RMS  Spot  Size  (pm) 

Optimized  Wavelet  Fit 
RMS  Spot  Size  (pm) 

-30.7° 

154.8 

18.7 

15.2 

-15.35° 

141.5 

11.5 

4.4 

0° 

191.1 

7.7 

4.2 

15.35° 

141.5 

13.5 

5.6 

30.7° 

154.8 

17.5 

10.7 

Table  2:  Both  fits  optimized  with  30  rays  per  bundle.  Spot  sizes  calculated  from  300  traced  rays  per  bundle. 

Diffraction  limited  spot  size  is  1.68  pm. 


B.4.c.  Some  results  that  show  how  much  better/faster  this  approach  can  be. 

The  results  thus  far  demonstrate  the  potential  of  basis  functions  as  tool  for  efficient  optimization  of  a  GRIN  lens, 
but  only  on  a  test  case.  To  demonstrate  the  design  engine  in  the  intended  application  a  TO/QCT  design  for  relay 
imager  was  used  in  the  second  study.  The  relay  imager  was  based  off  of  the  Luneburg  lens,  and  TO/QCT  was  used 


31 

Approved  for  public  release:  distribution  is  unlimited. 


to  make  the  input  side  concave  and  the  output  side  flat.  The  relay  imager  is  intended  to  be  used  a  monocentric 
camera.  The  initial  TO/QCT  design  for  the  imager  can  be  seen  in  Fig.  3. 
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Figure  3 


The  relay  imager  was  optimized  using  both  the  DCT  and  wavelet  bases  using  an  increasing  number  of  coefficients. 
The  largest  coefficients  were  used  for  this  simulation  given  that  the  initial  design  featured  sufficient  variation  in 
both  the  transverse  and  axial  directions.  The  behavior  of  this  lens  with  respect  to  the  number  of  optimization 
coefficients  was  far  less  predictable  than  the  quadratic  GRIN,  as  shown  in  Fig.  4.  Again  the  wavelet  basis  reached 
the  best  result  at  100  coefficients,  but  there  is  evidence  to  indicate  that  the  solutions  do  not  get  monotonically 
worse  with  increased  number  of  coefficients. 


Mo  no  centric  Relay 


Figure  4 

The  best  overall  design  is  shown  in  Fig.  5,  and  featured  spot  sizes  of  7.0,  8.9,  9.9,  7.8  and  3.4  microns  across  the 
field  of  view.  The  final  figure  of  merit  was  an  order  of  magnitude  improvement  over  the  initial  design. 
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Figure  5 

Due  to  the  obvious  dependence  on  initial  conditions  and  lens  geometry,  it  is  difficult  to  make  conclusive  and  far 
reaching  determinations  from  this  initial  study.  However,  it  was  clearly  demonstrated  that  the  basis  function 
approach  to  optimization  has  the  ability  significantly  improve  the  performance  of  the  TO/QCT  GRIN  lens  designs. 
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